The two-loop self-energy: diagrams in the coordinate-momentum representation 
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The paper reports a technique of evaluation of Feynman diagrams in the mixed coordinate- 
momentum representation. The technique is employed for a recalculation of the two-loop self-energy 
correction for the ground state of hydrogen-like ions with the nuclear charge numbers Z = 10 — 30. 
The numerical accuracy is considerably improved as compared to the previous calculations. The 
higher-order (in Za) remainder function is inferred from the numerical results and extrapolated 
towards Z — and 1. The extrapolated value for hydrogen is consistent (but still not in perfect 
agreement) with the analytical result obtained within the perturbative approach. 
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I. INTRODUCTION 

Investigations of the Lamb shift in atomic systems pro- 
vide one of the most stringent tests of quantum electro- 
dynamics (QED). They are also used for the determi- 
nation of fundamental physical constants [1]. The main 
factors limiting the present theoretical understanding of 
the Lamb shift are the binding two-loop QED effects and, 
first of all, the two-loop self-energy correction. 

Theoretical investigations of QED effects in light atoms 
traditionally rely on the perturbative expansion in the 
binding-strength parameter Za {Z is the nuclear charge 
and a is the fine structure constant). The state of the 
art of such calculations is the evaluation of the dominant 
part of the two-loop correction to order ma 2 (Za) 6 [2- 
5]. The main problem of the .Zen-expansion approach is 
the difficulty of estimation of uncalculated higher-order 
effects. In the case of the two-loop self-energy correction, 
the higher-order binding effects are above the experimen- 
tal error both for the light systems (particularly, for the 
hydrogen atom [6]) and for the heavy ions [7-10]. 

In the present investigation, we use the all-order ap- 
proach which is nonperturbative in the parameter Za. 
The nonperturbative calculations started with the pio- 
neering works of Wichmann and Kroll [11] and P. J. Mohr 
[12, 13]. For heavy ions, the all-order approach is the only 
alternative as the parameter of Za is of order of unity. 
For light systems, this method is complementary to the 
Za-expansion approach and can provide results for the 
high-order remainder beyond the known Za-expansion 
terms. 

The all-order calculation of the two-loop self-energy 
correction depicted on Fig. 1 was a long and difficult 
project accomplished in a series of papers [14-19]. The 
numerical results obtained in these studies agreed well 
with the first terms of the Za expansion calculated within 
the perturbative approach. However, a significant dis- 
agreement was reported [18] for the contribution to order 
ma 2 (Za) 6 (the so-called B^o coefficient). A reliable de- 
termination of this contribution from the all-order results 
requires a high numerical accuracy to be achieved for the 
low values of Z, which is a challenging task. In our recent 



investigation [20] , we briefly reported a new calculational 
technique for the evaluation of Feynman diagrams in the 
mixed coordinate-momentum representation. This tech- 
nique significantly improved the numerical accuracy of 
the two-loop self-energy calculation and to a large extent 
removed the disagreement with the analytical approach. 
In the present paper, we present a detailed description of 
this calculational technique. 

The relativistic units (m — h — c — 1) and the Heav- 



iside charge units (a 
throughout the paper. 



e /47r, e < 0) will be used 



II. TWO-LOOP SELF-ENERGY 

The Feynman diagrams representing the two-loop self- 
energy correction are shown in Fig. 1. The contribution 
of the first diagram [Fig. 1(a)] is conveniently divided 
into two parts, the irreducible and the reducible one. The 
reducible part is induced by the virtual states with the 
energy e„ = e a in the middle electron propagator (e Q is 
the energy of the reference state), and the irreducible part 
is the remainder. The irreducible part (often referred to 
as the loop-after-loop correction) can be interpreted as a 
second-order perturbation induced by the one-loop self- 
energy operator. The corresponding expression reads 
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where S(e a ) = £(e a ) — Sm, S(e) is the one-loop self- 
energy operator [17], 5m is the corresponding mass coun- 
terterm, and G rod is the reduced Dirac-Coulomb Green 
function. The irreducible part is finite and can be calcu- 
lated separately by generalizing various methods devel- 
oped for the one-loop self-energy. 
The reducible part is given by 



A£ red = A£sE (7°^S(e) 
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where AEse — (7°S(e a )) is the one-loop self-energy cor- 
rection. 
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FIG. 1: Feynman diagrams representing the two- loop self- 
energy correction in the external binding field. The double 
line represents the electron propagating in the field of the 
nucleus. 



The contribution induced by the diagram in Fig. 1(b) 
will be referred to as the overlapping term. It is given by 

AEo = 2ia / duj\ / dx\ . . . dx± D(u>i, £13) ipl(xi) 
J —00 J 

x G(e a - wi) 7°A A '(e a - wi,e a ) ip a (x4) - Sm , 

(3) 

where Smo is the mass counterterm, D(u>,xi 2 ) is the ra- 
dial part of the photon propagator in the Feynman gauge, 



D(lj, X12) = 



exp 



(i Vw 2 + iOxi 2 ) 



4irx 
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(4) 



G(e) is the Dirac-Coulomb Green function defined by 
G(e) = [e— H(l— iO)]" 1 , with H being the Dirac Coulomb 
Hamiltonian, and X12 = \x\ — x 2 \. The vertex function 
A M is defined as 

/oo 
-00 

x G(e a - u x - lj 2 ) a" G(e a - oj 2 ) a v . (5) 

The contribution induced by the diagram in Fig. 1(c) 
will be referred to as the nested term. It reads 

AEn — 2ia / duj\ I dxi . . . dx± D(u>i, Xu) ipl(xi) 
J — 00 J 

x a^G(s a - bJi) 7°S(e a - oji) G(e a - Wi) ip a (x 4 ) 



(6) 



where Stun denotes the mass counterterm. 

The general analysis [17] shows that the sum of the 
reducible, the overlapping, and the nested terms is finite. 
However, the individual contributions are divergent both 
in the ultraviolet and the infrared regions of virtual pho- 
ton energies. In order to make all contributions explicitly 
finite and suitable for a numerical evaluation, a careful 
rearrangement of individual parts is required. This re- 
arrangement is discussed in detail in Rcf. [17] and will 
not be repeated here. The general idea is that the bound 
electron propagators in the loops are expanded in terms 
of the interaction with the binding field and the result- 
ing contributions are grouped together into three large 
classes: (i) the part calculated in the coordinate space, 
conventionally termed as the M term and denoted by 
AEm, (ii) the part calculated in the momentum space 
(the F term AEp), and (iii) the part calculated in the 
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FIG. 2: The P term. The single line represents the free elec- 
tron propagator. The dashed line with a cross indicates the 
interaction with the Coulomb field of the nucleus. 



mixed coordinate-momentum representation (the P term 
AE P ), 



AE Icd + AE + AE N = AE M + AE F + AE P . 



(7) 



All the three terms can be made explicitly finite and cal- 
culated separately. The calculational technique is com- 
pletely different for each term. In the present investi- 
gation, we concentrate on the P term, as the scheme of 
evaluation of the other two was described in detail in 
Rcf. [17] and has not been changed significantly since 
that work. 



III. P TERM: BASIC FORMULAS 

The Feynman diagrams contributing to the P term are 
shown in Fig. 2. They arise from the diagrams in Fig. 1 
when the bound-electron propagators in the loops are ex- 
panded in terms of the interaction with the binding field. 
The characteristic feature of the diagrams on Fig. 2 is 
that the ultraviolet divergences in them originate from 
the one-loop subgraphs only. The divergent subgraphs 
are covariantly regularized and calculated in the momen- 
tum space, whereas the remaining part of the diagrams 
does not need any regularization and is evaluated in the 
coordinate space. 

It should be mentioned that the necessity of calcu- 
lation of Feynman diagrams in the mixed coordinate- 
momentum representation is a distinctive feature of the 
two-loop effects treated to all orders in the parameter 
Za. Because of this, the P term has no analog nei- 
ther in the one-loop calculations nor, to the best of our 
knowledge, in any other previous calculations. (All other 
two-loop effects evaluated to all orders in Za were ef- 
fectively reduced to one-loop contributions, see Ref. [21] 
and references therein.) For the first time the P term 
was calculated in Rcf. [15] with help of a finite basis set 
representation of the spectrum of the Dirac equation. In 
the present investigation we report a different technique 
based on the analytical representation of the Green func- 
tion in terms of the Whittakcr functions. 
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As shown in Fig. 2, the P term is represented by a sum 
of five terms, 

AE P = AEp. a + AE P , b + AE P , C + AE P . d + AE P . e , 

(8) 

each of which refers to the corresponding diagram. In or- 
der to make the individual terms finite, we assume that 
the one-loop subgraphs are represented by the renormal- 
ized operators and that the infrared reference-state singu- 
larities are removed by the minimal subtractions. It can 
be explicitly checked that such definition of the P term 
is equivalent to the definition of Ref. [17], so that the 
present numerical results are directly comparable with 
those of the previous work. 

The contribution of the diagram in Fig. 2(a) is given 

by 

AEp, a = 2iaJ dut J J dxidx 2 D(u, x\ 2 ) 

x Vl(a;i)a M G v (E,x 1 ,p)S 1 (E,p)G v (E,p,x 2 ) 

- G ( « } (E, Xl ,p) Si (e a ,p) G^ (E, p, x 2 )] a^a{x 2 ) ■ 

(9) 



where E = e„ 



1 



S1(S,P) = Q 

7 u e — 7 • p — m 



7 U £ — 7 • p — m 

( 10 ) 

is the renormalized free self-energy operator (see 
Ref. [17] for its definition and explicit representations), 
and the function Gy is the Fourier transform of the prod- 
uct of the Green function G and the Coulomb potential 

V c (x) = -Za/x, 



G v (s,xi 



dx 2 e ip - X2 G{e,x 1 ,x 2 )V c {x 2 ), (11) 



G v (e,p,x 2 )= j dx 1 e- l P xl V c {x 1 )G{e,x l ,x 2 ). (12) 

The second term in the brackets of Eq. (9) removes the 
reference-state singularity present in the first term. The 
definition of the function Gy ' is obtained from Eqs. (11) 
and (12) by the substitution G -> G< a \ where G^ is the 
reference-state part of the electron propagator defined by 

^)( e , 81 , 8a) = ^ ^gjj^M . (13) 



Here, tp a i denotes the virtual electron state of the same 
energy and of the same parity as the reference state and 
\x a i is its momentum projection. 

The contribution of the diagram in Fig. 2(b) is given 



by (with the combinatorial factor of 2) 

doj J ^2^)3 J dx\dx 2 D(u>,xi 2 ) 
x ^(xi) a M [g v (E, x u p)- G { °\E, xi, p) 
x S 2 (E,p) G (0 He,p, x 2 ) a" M*2) , (14) 

where 



4 0) (£,P), (15) 



y'e — 7 • p — m 



G^ is the free Dirac Green function, and the function 
Gy^ is defined by Eqs. (11) and (12) after the substitu- 
tion G->G(°). 

The contribution of the diagram in Fig. 2(c) is 

AE P , C =2iaJ_ oo d U /(|y3 (Sj3 j ***** 
x D{uj,x 12 )V c {q) ^iiix^a^ 
x G v (E,x 1 ,p 1 )g 1 (E,p 1 ,p 2 )G v (E,p 2 ,x 2 ) 

-G ( v l) (E,x 1 ,p 1 )G 1 (e a ,p 1 ,p 2 )G { * ) (E,p 2 ,x 2 ) 
xa^ a (x 2 ), (16) 



where Vc(q) = — AnZa/q 2 is the Coulomb potential in 
momentum space, q = Pi — p 2 , 



1 



7 u e — 7 • p± — m 
x r^(e,pi;£,p 2 ) — 



7 u e — 7 • p 2 — m 



(17) 



and T° R is the time component of the renormalized free 
vertex operator T^, (its explicit representation can be 
found in Ref. [17]). 

The contribution of the diagram in Fig. 2(d) is given 
by (with the combinatorial factor of 2) 



AEt 



4/n / duj I ,J } ] n ,y 2 ^ I dx\dx 2 



(2tt) 3 (2tt) 3 
x D(lj,x 12 ) V c (q)tpl(xi)a t _ t G v (E,Xi,pi) 
x g 2 (E,p u p 2 )G^(E,p 2 ,x 2 )a^i; a (x 2 ), (18) 



where 



1 



T° R (e, Pl ;s,p 2 ). (19) 



5 2 (e,pi,p 2 ) = — 

7 u e — 7 • pi — m 

The contribution of the diagram in Fig. 2(e) can be 
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written as (with the combinatorial factor of 2) 

dpi dp 2 



AE Pe = -Aia I dio 

J — oo 

exp(— iq ■ z) 
uj 2 — q 2 + iO 



j°E — 7 • px — m 



(2tt) 3 (2tt) 3 ' <IZ 

(0) ( 

T R (E,p 1 ;e a ,p 2 )'ip a (p2) ■ 



G v (E 7 z 7 p 1 )-G^(E,z,p 1 



(20) 



IV. CALCULATION 



The main problem of the evaluation of the P term is 
connected with the fact that the regularization of ultravi- 
olet divergences in the one-loop subgraphs is carried out 
in the momentum space, while the bound-electron propa- 
gators are most easily evaluated in the coordinate space. 
As a result, the expressions listed in the previous section 
contain the Fourier transform of the product of the Dirac 
Coulomb Green functions G with the Coulomb potential 
Vc over one of the radial arguments, see Eqs. (11) and 
(12). Since we were not able to find a satisfactory analyt- 
ical representation for such an object, the only possible 
way was to evaluate the Fourier tansform integral nu- 
merically. This way entails numerous numerical integra- 
tions of rapidly oscillating functions, which may make the 
computations prohibitively expensive, unless great care 
is taken in choosing the optimal calculational approach. 

In the previous investigations [15-17], the P term was 
calculated with help of a finite basis set for the Dirac 
equation. The advantage of the basis-set methods is that 
they represent the Dirac Green function as a continuous 
function of the radial arguments (for any hnite size of the 
basis), whereas the exact Green function is discontinous 
when the two radial arguments are equal. The usage 
of a finite basis set allows one to perform the Fourier 
transform of the Green function over one radial variable 
independently of the other. However, the convergence 
with respect to the size of the basis appears to be the 
limiting factor for the accuracy of the calculations. In 
the present investigation, we set up a different calcula- 
tional approach based on the analytical representation of 
the Green function. Technical details of evaluation of the 
Dirac Coulomb Green function in the mixed coordinate- 
momentum representation are described in Appendix A. 
The corresponding formulas for the free Dirac Green 
function are summarized in Appendix B. 



A. Nested diagrams 

In this subsection, we address the diagrams in 
Figs. 2(a)-(d) and outline the major steps required to 
make the basic formulas suitable for a numerical evalua- 
tion. 



The integrations over the angular variables [x\, £2, 
and p in Eqs. (9) and (14) and x\, x 2 , pi, and p 2 in 
Eqs. (16) and (18), where x — x/\x\] are relatively sim- 
ple. We employ the fact that the matrix elements of the 
operators S1.2 and Qi i2 with the Dirac wave functions 
are (i) diagonal with respect to the angular momentum 
quantum number k and the momentum projection /j, and 
(ii) do not depend on /i: 

(Kfj,\Si\K' fj,') = 6 KtK > 5^ {n\Si\K) , (21) 
MGi Vbl/sV) = <W<W (K\QiV C \K) . (22) 

As a result, the integrations over p in Eqs. (9) and (14) 
and pi and p 2 in Eqs. (16) and (18) are exactly the same 
as for the zero- and one-potential parts of the one-loop 
self-energy correction, see Ref. [22] for details. The inte- 
grations over Xi and x 2 are the same as for the many- 
potential part of the one-loop self-energy correction. 

As illustrated in Ref. [22] , the intcgations over p\ and 
P2 in Eqs. (16) and (18) can be easily reduced to a sin- 
gle integral over £ = Pi ■ P2, which needs to be evalu- 
ated numerically. The calculation is complicated by the 
presence of an integrable Coulomb singularity (~ l/<? 2 , 
q = \pi — p 2 \). This singularity is removed in two steps. 
First, the change of the integration variable £ — > q weak- 
ens it to ~ l/q. The remaining singularity is removed by 
subtraction of the vertex function with the equal argu- 
ments, 

t° r (jpi,P2) -> T° R {p uP 2) - \ [rSj(pi,pi) + r° R (P2,P2)] . 

The vertex operator with two equal arguments is related 
to the free self-energy operator by the Ward identity: 



« = -^'(p)- 



(24) 



In the subtracted terms, the Coulomb singularity is eas- 
ily integrated out by using identities obtained from the 
definition of the Dirac Coulomb Green function, as, e.g., 



dp2 



^(</)-7n 



1 



■G v {E,p 2 ,x 2 ) 



(2ir) 3 ' ~ v "* / 7 E — 7 • P2 — m 

G v {E, Pl ,X2)-Gf{E,p u X2). (25) 



Finally, we change the contour of the integration over 
the virtual photon energy ui in Eqs. (9), (14), (16), and 
(18) from (—00,00) to a new contour Cl#, whose main 
part is parallel to the imaginary axis. The contour Clh 
consists of the low-energy Cl and the high-energy Ch 
parts. The low-energy part extends over (A — iO, —iO) 
on the lower bank of the cut of the photon propagator 
and over (iO, A + iO) on the upper bank of the cut, with 
the parameter A fixed by A = Zae a - The high-energy 
part consists of two intervals, (A + iO, A + ioo) and (A — 
iO, A — ioo). The contour Clh differs from the one used 
by P. J. Mohr [12] only by the choice of the breaking 
point A (the value A = e a was employed in that work). 



■5 



In our previous calculations we used the contour that 
extended along the imaginary axis (which corresponds 
to the choice of A = 0). The contour Clh is more con- 
venient for the numerical evaluation. Firstly, there is 
no pole contributions originating from the reference-state 
part of the electron propagators. Secondly and more im- 
portantly, the photon propagator on the low-energy part 
of the contour involves sin(wa;i2), which suppresses small 
denominators due to the virtual bound states and leads 
to a smooth behaviour of the integrand for small uj. 

In the present investigation, we are concerned with the 
reference state being the ground state only. In this case, 
no pole contributions appears for the contour Clh- For 
the excited reference states, however, there are single and 
double poles on the low-energy part of the contour, which 
arise from the intermediate states bounded more deeply 
than the reference state. These poles require a separate 
treatment or a deformation of the integration contour 
into the complex plane. 



B. Overlapping diagram 

In this subsection we address the overlapping diagram 
shown in Fig. 2(e), whose expression is given by Eq. (20). 
The angular integration in this expression is rather in- 
volved and will be considered in detail. 

In order to perform the integration over z, we expand 
the exponent into the spherical waves, 

e-** = 4nJ2^ L JL(qz)Y LM (q)Yl M (z) , (26) 

LM 

where jl is the spherical Bessel function and Ylm is the 
spherical harmonics. The time component (p = 0) of the 
z integration is immediately evaluated in terms of the 
basic integrals of the form 

J dz X i b ^(z)Y LM (z) XKa , a (z) = S |V (K b \\ C< L >||« > , 

(27) 

where Xkh(z) are the Dirac spin- angular spinors [23], 
is the spherical tensor with components Cff {$) = 
y/4n/(2L + 1) Ilm(^), (II ' ' ' II) denotes the reduced ma- 
trix element, and 



ba _ 
S LM — 



(-1) 



47T 



C 



LM 

3b ,3a- Ha 



(28) 



with Cj^ ni j 2rri2 being the Clebsch-Gordan coefficient. 

The vector components (fi = 1, 2, 3) of the z inte- 
gration are calculated after expanding the integrand in 
terms of the vector spherical harmonics Yjlm [24, 25], 

xt bfib (z) <rX Ka na(z) = S JM SjL(K b ,K a )Yj LM (z) , 



JLM 



where er is a vector incorporating Pauli matrices. The 
coefficients Sjl are given by 



Sjj+i(n a , n b ) 



S.Jj(n a ,K h ) 



SjJ-l{K a , K b ) = 



J +1 I ^ K„ ' Hi, 



2J+ 1 V J+l 

X (- Kb \\C^\\K a ), 



(30) 



K a Kb 

y/J(J+ 1 ) 



J 



(K b \\C^\\ Ka ), (31) 



1 



2J+ 1 V J 

x (-K b \\C^\\ Ka ). 



(32) 



For J = 0, the only nonvanishing coefficient is Soi- 

We now turn to the evaluation of the integrals over pi 
and p 2 in Eq. (20). The aim is to integrate out all angu- 
lar variables except £ = pi ■ p 2 . To this end, we exam- 
ine the angular structures encountered in the integrand. 
The time component of the vertex operator sandwiched 
between the Dirac wave functions involves two indepen- 
dent angular structures, which are identified by 



^(Pi)-o 



j°E — 7 • pi — m 



^j\(E,Pi]e a ,p 2 )ip a (p 2 ) 



= ^ ^ la { [9n Fig + fn Flf] xt n „ n (Pi) X*aHa (&) 
+ [gn^g + fnftf] X- Knftn (Pl) X-K aMa (p 2 )| , (33) 

where g n = g n (pi) and /„ = f n (pi) are the upper and 
the lower radial components of the Dirac wave function 
ip n and Ti are scalar functions Ti = J r i(E,s a ,pi,p2,q), 
with pi = \pi\, p 2 = \p 2 \, and q = \q\. The vector part 
of the vertex operator induces six angular structures, 



T R (E,p 1 ;e a ,p 2 )tp a (p2) 



(29) 



= {[9nK lg + f n lZ lf ] XL n ^ n {Pl)^X-KaHa(P2) 

+ [9n K 2g + f n TZ 2f ] X-n n a n (Pi) <* Xn^ a (p 2 ) 
+ [gn n 3g + f n Tl 3 f] pi xi n a n (Pl)X^a (p2 ) 
+ [gn K Ag + f n TZ 4 f] p 2 xi nlin {Pl)X^ a n a (j5 2 ) 
+ [gn K 5g + f n H 5 f] Pl X- Kn a n {Pl)X-K a u a (p 2 ) 

+ [g n n 6g + f n 1lef] P2X- K „ M „(pi)x- KaM a(p2)} , (34) 

where IZi = lZi(E,e a ,Pi,P2,q)- The functions Ti and 
IZi can be straightforwardly obtained from formulas in 
Appendix A of Ref. [26]. 

Using Eqs. (33) and (34), it is possible to parameter- 
ize the angular structure of the integrand of Eq. (20) by 
four basic angular factors t KltK2 , s% 1>K2 , s p K \ iK2 , and S P 2 i K2 
defined as 

t Kn , K a(J) = J2 s Tm xi n „ n {Pi)Y JM {q) XwM , 

ti n M 

(35) 
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(36) 

«?.,».(^) = J2 s 7Mxi n ^(Pi)P t -YjLM{q)x^ a {P2)- 

fl n M 

(37) 

By an explicit evaluation with help of formulas from 
Ref. [25] , one can show that the above angular factors are 
the functions of pi, p 2 , and q only (or, in other words, 
that they depend on the angular variables only through 
£). This statement allows us to integrate out all angles 
in Eq. (20) except £. The calculation of angular factors 



is described in Appendix C. 

For a numerical evaluation of Eq. (20), we need to de- 
form the contour of the co integration. In this case (in 
contrast to the nested contributions), we find it conve- 
nient just to rotate the integration contour to the imagi- 
nary axis, lj — > iuj. This leads to appearance of the pole 
contribution. So, 

AE P . e = AE^ e + AE p P f , (38) 

where AEp^ e is the contribution of the integral along the 
imaginary axis and AE p ° e e is the pole contribution. 

The final result after the angular integrations and the 
rotation of the contour is 



2 roo poo p 

AE l ™ e = -r 5ft y / cLu / dpi dq / 
^ „ Jo Jo J\t 



Pl+q 

dpi dq I dp 2 

\Pi-Q\ 



qpiPi 

-w 2 - q 2 J 



dz z 



x |^(-l) fa i7M (n n \\C^\K) [{9aG^ n +f a G 2 JjF g + (g a G™ n +f a GfjFf 
- Y,(-l) k2 3L(qz) ([«?„ d$ Hn S. JL (n a ,-n n ) - f a G v \ n S JL (- Ka ,K„)] K g 

JL 

+ [da Gy Kn S JL(Ka, -K n ) - f a Gy^ S JL (-K a , K n )] TZf^j j , 



(39) 



where h = ( J - l n + l a )/2, k 2 = (L — l n + l a — l)/2, 
l n = |«n + 1/ 2 I - V 2 > .9a = ffa(^) and / a = f a (z) are the 
radial components of the reference-state wave function, 
and Gy stand for the difference of the radial compo- 
nents of the Coulomb Green function (times Vc) and the 
free Green function (times Vc), 

Gy Kn = G VKn (e a - ioj, z,pi) - Gy>^ (s a - ioj, z,p x ) . 
The angular functions in Eq. (39) are defined by 

T g = T lg t Kn , Ka (J) + T 2g *-„„,_„„ (J) , (40) 

K g =1Z lg s° n! _ Ka (JL)+1l 2g s°_ K ^ Ka (JL) 
+ PiK 3g sZ, Ka (JL) + p 2 TZ Ag 8* tKa (JL) 
+ p 1 n 5g s p } Knt _ Ka (JL) + P2 Tl 6g s p _\ nt _ Ka (JL), (41) 

and the same for the J 7 / and IZf functions. 

The pole contribution A£P° lc is obtained from Eq. (39) 
by the following substitution (valid for a being the ground 
state), 

-+ -\ S KnKa 5{u>) cPi(z) <jP Va ( Pl ) , (42) 



where 4>l(z) = g a (z), <\> 2 a (z) = f a (z), and <p Va (p) is the 
Fourier transform of the product (x) Vc (x) . 

A useful check of the angular-momentum algebra con- 
sists in making the substitution Tj i (E,pi;e a ,p 2 ) — > 7 P 
in Eq. (20). The result is an one- loop self-energy con- 
tribution which can be calculated independently in the 
coordinate representation. 



C. Numerical evaluation 

We start our discussion of the numerical evaluation of 
the P term with AEp :a and AEp^ given by Eqs. (9) and 
(14), respectively. These are the two simplest contribu- 
tions. After carrying out integrations over the angular 
variables as described in Sec. IV A, five integrations re- 
main to be calculated numerically, namely those over w, 
p, x\, and x 2 and the Bessel transform integral implicitly 
present in the Green function. 

The radial integrations over x\ and x 2 have to be or- 
ganized in such a way as to avoid unnecessary recalcu- 
lation of the Bessel transform integrals, as discussed in 
Appendix A. To this end, we set up a radial grid {xij t k} 
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as follows. The first-level elements Xj ; o,o are given by 

l_f? 

a;»,o,o = Po 2 1 , (43) 

z 

where p is a parameter adjusted empirically, the vari- 
able ti is uniformly distributed over the interval (t m in, 1), 
and a small value of i m ; n > cuts off the radial inte- 
grations at large distances. The second-level elements 
Xijfi represent the Gauss-Legendre quadratures on the 
interval (xifi,o, ^<+i,o,o)- The third-level elements Xij^k 
represent the Gauss-Legendre quadratures on the inter- 
val (xijfi,Xi,j+i,o)- In the result, we obtain an ordered 
three-level radial grid. To perform the radial integrations 
over X\ and X2, it is sufficient to know the integrand on 
this grid only. 

The general scheme of the evaluation of AEp^ a looks as 
follows. For fixed values of n, lo, and p, we set up the ra- 
dial grid {xij t k}- Ori this radial grid, we store the compo- 
nents of the Dirac Green function (jP K {E, x) and (f>^(E, x) 
[see Eqs. (A3) and (A4)], the Bessel transform functions 
ip° K (E,p;x) and ijj™(E,p;x) [see Eqs. (A9) and (A10)], 
and the other functions required for the evaluation of the 
integrand (the radial part of the photon propagator, the 
reference-state wave function, etc.). After that, the ra- 
dial integrations are performed simply by summing up 
the stored numerical values. Next, we perform the in- 
tegration over p, then the one over w, and finally, the 
summation over n. 

The most expensive part of the calculation is the eval- 
uation of the Bessel transform integrals. In order to con- 
trol the accuracy of numerical integrations, one needs 
an efficient procedure for calculating the transforms for 
various momenta, including ones as large as 10 6 . In 
our calculations, we used the Gauss-Legendre integra- 
tion quadratures in the region where the argument of 
the Bessel function is of about unity or smaller. Outside 
this region, the spherical Bessel function was expressed 
as a combination of the sine and cosine functions. The 
numerical evaluation of the sine and cosine transform was 
performed with help of routines of the NAG Fortran li- 
brary. 

The scheme described above works well for the AEp^ a 
contribution but turns out to be not sufficiently effec- 
tive for AEp.t,, leading to a slow convergence of the 
radial integrations with respect to the number of inte- 
gration points. This is because the free Green function 
G^°\e,x,p) contains a Bessel function [see Eq. (B5)], 
which oscillates rapidly in the high-momenta region. 
This problem was solved by observing that the integral 
over X2 in Eq. (14) has a structure similar to tl> K (E,p; xi), 
i.e., it is essentially a Bessel transform over the intervals 
(0,21) and (21,00). We thus perform the integration 
over X2 in Eq. (14) by means of the same approach as 
used in the evaluation of the Bessel transform functions 
ip K (E,p; x\). This approach improves the convergence of 
the radial integrals drastically. 

The evaluation of the two remaining nested contribu- 
tions, AEp c and AEp^, was performed in the full anal- 



ogy with the discussed above. However, it turned out to 
be much more time consuming due to a larger number 
of integrations. Indeed, in place of an integration over 
p in AEp. a and AEp^, there are now four integrations 
(those over p\ , p 2 , and q and the Feynman-parameter in- 
tegration implicitly present in the vertex operator). For- 
tunately, the integrations over q and over the Feynman 
parameter can be carried out independently of the inte- 
grations over x\ and x 2 and thus do not significantly in- 
fluence the total calculational time. The two integrations 
over pi and P2, however, lead to a considerable increase 
of the computational expense (this being about a week 
of the processor time for each value of Z). 

For the overlapping contribution AEp^ e , there are 
seven integrations to be performed numerically Five of 
them are explicitly written in Eq. (39), one over the Feyn- 
man parameter is implicitly present in the vertex opera- 
tor, and the last one is the Bessel transform integral in 
the Green function. The number of nested integrations 
can be reduced by observing that the integrations over p 2 
and over the Feynman parameter can be carried out in- 
dependently of the integration over z. The calculation is 
complicated by the fact that the integral over z contains 
a Bessel function, thus being essentially a Bessel trans- 
form. In order to get a stable result for the z integration 
in the region of large values of qz, we had to interpo- 
late the part of the integrand that multiplies the Bessel 
function and evaluate the Bessel transform analytically. 



V. RESULTS AND DISCUSSION 

The results of our calculation of the P term for the 
ground state of hydrogen-like ions with Z > 10 are listed 
in Table I. The individual contributions are presented 
in a way that allows a detailed comparison with the pre- 
vious calculations. Specifically, the first four columns of 
Table I are directly comparable to the four columns of Ta- 
ble 2 in Ref. [17]. The previous results listed in the fifth 
column of Table I were obtained in Ref. [18] for Z < 60 
and in Ref. [17] for the other Z. The agreement with 
the previous calculations is very good in most cases, but 
the present numerical accuracy is significantly higher. A 
small deviation in the high-Z region is probably due to a 
difference in the treatment of the nucleus. (In the present 
work, the point nuclear model is used, whereas the previ- 
ous investigations [17, 18] were conducted with a partial 
inclusion of the finite nuclear size effect in the P term.) 

The uncertainty of the present results is mainly due 
to the termination of the partial-wave expansion. In our 
calculation, we included typically 20-30 partial waves and 
estimated the omitted tail by fitting the data obtained as 
a function of the cutoff parameter. Perspectives for im- 
proving the present accuracy further (which is required 
if one is to perform a calculation for lower values of Z) 
seem to be questionable. The main problem is that the 
high partial waves become increasingly difficult to con- 
trol numerically. At the same time, the extrapolation of 



8 



TABLE I: The P term for the ground state of hydrogen-like ions, in units of AE/[ma 2 (Za) 4 /iv 2 ]. 



Zj 


Figs. (a,b) 


Figs. (c,d) 


£ig. (e) 




previous [if, loj 


10 


-855.4289 (20) 


1265.4550(50) 


-1131.3372 (22) 


-721.3111 (58) 


-721.34(12) 


12 


-486.0740 (20) 


744.6950 (50) 


-697.6864 (17) 


-439.0654 (56) 




15 


-239.2533 (15) 


384.3590 (30) 


-380.3170(12) 


-235.2113 (36) 


-235.205 (70) 


17 


-159.4480(15) 


263.5625 (20) 


-268.3945 (12) 


-164.2800(28) 




20 


-93.3597(7) 


160.3551 (20) 


-169.0247(10) 


-102.0293(23) 


-102.026 (55) 


25 


-44.1892(7) 


79.9950 (15) 


-87.7882(10) 


-51.9824(19) 




30 


-23.8155 (7) 


44.8102 (27) 


-50.4083 (10) 


-29.4135 (30) 


-29.410 (25) 


40 


-9.0138 (4) 


17.6061 (2) 


-20.1713 (6) 


-11.5790(7) 


-11.575 (30) 


50 


-4.3391 (5) 


8.4020 (4) 


-9.5506 (4) 


-5.4877 (7) 


-5.488 (26) 


60 


-2.4455 (3) 


4.5451 (4) 


-5.0660 (2) 


-2.9664 (5) 


-2.970(18) 


70 


-1.5203 (2) 


2.6716(1) 


-2.9354(2) 


-1.7841 (3) 


-1.757(25) 


83 


-0.8655 (2) 


1.4268(1) 


-1.6307(1) 


-1.0693 (2) 


-1.057(13) 


92 


-0.5545(1) 


0.9091 (1) 


-1.1902(1) 


-0.8356 (2) 


-0.812(10) 


100 


-0.2990 (2) 


0.5426(1) 


-0.9792 (4) 


-0.7356 (5) 


-0.723 (7) 



the partial wave expansion requires an accurate represen- 
tation of the individual partial-wave expansion terms. 

The main motivation of the present investigation was 
to improve the numerical accuracy of the total two-loop 
self-energy correction in the region of medium values of 
Z, in order to get a more reliable extrapolation towards 
Z = 0. To this end, a new approach was developed for 
the evaluation of the P term, as described above. Besides 
that, the other parts of the two-loop self-energy correc- 
tion were reevaluated to a higher accuracy. This was ac- 
complished by the methods described in Ref. [17], with 
the increased number of partial waves included and with 
denser integration grids. The results were first presented 
in Ref. [20]. 

Table II summarizes our results for the total two-loop 
self-energy correction for the ground state of hydrogen- 
like ions with the nuclear charge Z = 10 — 30. We observe 
that the calculational errors of the P term do not influ- 
ence significantly the errors of the total results, as the 
main uncertainty is now delivered by the M term. This 
uncertainty originates both from the dependence of the 
results on the number of integration points and from the 
termination of the partial- wave expansions. Since there 
are two independent partial-wave expansion parameters 
in the M term (see Ref. [17] for details), the number of 
expansion terms grows drastically as the cutoff parame- 
ter is increased. Because of this, significant extension of 
the partial-wave summations looks prohibitively expen- 
sive at present. 

The two-loop self-energy correction for the ground 
state of hydrogen-like atoms can be conveniently repre- 
sented in the following form, separating out the known 
terms of the Za expansion, 



AE = m g) 2 (Za) 4 [s 4 o + (Za) G 50 (Z)] , (44) 



and 



(Zar 2 ]B t 



63 



G 50 (Z) =B 50 + (Za){ln 3 [(l 

+ \n\Za)- 2 ] B 62 + \n[(Za)- 2 ] B el + G 60 (Z)} , 

(45) 



where B^ are the expansion coefficients with the first 
index corresponding to the power of Za and the sec- 
ond index, to the power of logarithm, and Gij(Z) are 
the functions incorporating the corresponding Bij and 
all higher orders in Za, Gij(Z) = B^ + Za (■■■)■ The 
results for the expansion coefficients (see Refs. [1-5] 
and references therein) are: B40 = 1.409244, B50 = 
-24.2668(31), B 63 = -8/ 27, B 62 = 16/ 27 - (16/9) In 2, 
B 61 = 48.388913, and B 60 = -61.6(9.2). 

The functions Gij(Z) inferred from our numerical data 
are plotted in Fig. 3. The visual agreement of our results 
with the analytical values of the expansion coefficients is 
very good for B40 and -B50, but not exactly satisfactory 
for B eo . In order to produce a clearer statement, we need 
to extrapolate our data towards low values of Z. For 
this we use a variant of the procedure first employed in 
Ref. [27]. The extrapolation towards the required value 
of Z = Zo (=0 and 1 in our case) is performed in two 
steps. First, we apply an (exact) linear fit to each pair 
of two consecutive points from our data set and store the 
resulting values at Z = Z - Second, we perform a global 
parabolic least-squares fit to the set of data obtained on 
the first step and take the fitted value at Z = Zq as a 
final result. Similar procedure applied to the function 
G50 reproduces the analytical result for the coefficient 
B50 with the accuracy of about 1%. For comparison, a 
global polynomial fit yields a result accurate within 5% 
only. 

When applied to the remainder function Gqq(Z), the 
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TABLE II: The two-loop self-energy correction for the ground state of hydrogen-like ions, in units of AE /[ma 2 (Za) 4 /it 2 ]. 



z 


LAL 


F term 


P term 


M term 


Total 


2005 results [18] 


10 


-0.358 


822.138 (5) 


-721.311(6) 


-100.297(35) 


0.172 (36) 


0.25(16) 


12 


-0.417 


519.603 (2) 


-439.065 (6) 


-80.117(38) 


0.004 (38) 




15 


-0.495 


292.901 (2) 


-235.211 (4) 


-57.406(11) 


-0.212(12) 


-0.164(85) 


17 


-0.541 


211.052 (1) 


-164.280 (3) 


-46.567(9) 


-0.336(10) 




20 


-0.602 


136.909 (1) 


-102.029 (2) 


-34.780 (4) 


-0.501 (5) 


-0.481 (58) 


25 


-0.686 


74.501 (1) 


-51.982 (2) 


-22.560 (6) 


-0.728(6) 


30 


-0.756 


44.728 (1) 


-29.414 (3) 


-15.468 (3) 


-0.910(5) 


-0.903 (26) 




5 10 15 20 25 30 5 10 15 20 25 30 5 10 15 20 25 30 



Z Z Z 



FIG. 3: The two-loop self-energy correction. Gao(Z) = AE /[ma 2 (Za) 4 /tt 2 ], G 50 {Z) and G 60 {Z) are defined by Eqs. (44) and 
(45). The cross on the y-axis indicates the analytical results, G<io(0) = B40, Gso(0) = B50, and G6o(0) = Beo- 



extrapolation procedure described above gives 

G 60 (Z = 0) --84(15), (46) 
G eo (Z = l) =-86(15). (47) 

The extrapolated value for Z = 1 is higher than but 
marginally consistent with the 2005 result of —127(42) 
[18]. The shift of the central value is due to two reasons. 
First, the analytical result for the Bq\ coefficient was 
recently changed by 8Bqi = —1.4494 ... [5], thus pushing 
the remainder function higher up. Second, the improved 
numerical accuracy of the present calculation and the 
increased number of values of Z studied allowed us to 
identify the upward trend in the numerical data. 

The shift of the extrapolated values for Geo signifi- 
cantly reduced the disagreement with the analytical cal- 
culation [3] reported in Ref. [18]. The present value of 
G6o(0) = —84(15) is consistent (but still not in perfect 
agreement) with the analytical result of B 60 = —62(9) . 

To complete our analysis of the higher-order two- 
loop effects in hydrogen, we combine the result for the 
two-loop self-energy correction obtained in this work 
[Eq. (47)] with the corresponding contribution induced 
by the two-loop diagrams with the closed fermion loops 
reported in Ref. [21]. So, our estimate of the total two- 
loop (nonlogarithmic) contribution to order ma 2 (Za) 6 



for the ground state of hydrogen is 

G 6Q (Z = 1, total) = -86 (15)- 15 (2) = -101(15). (48) 

The numerical contribution of this effect is 
— 10.2(1.5) kHz, which is much larger than the er- 
ror of the experimental determination of the IS* — 2S 
transition frequency in hydrogen [6] (34 Hz) and com- 
parable with the experimental errors for the 2S — 12D 
transitions [28] (7 kHz). 

To conclude, in the present investigation, we described 
in detail the technique of evaluation of Feynman dia- 
grams in the mixed coordinate-momentum representa- 
tion based on the analytical representation of the bound 
electron propagators in terms of the Whittaker functions. 
This technique allowed us to significantly improve the ac- 
curacy of the numerical evaluation of the part of the two- 
loop self-energy correction conventionally termed as the 
P term. The all-order (in the parameter Za) results are 
reported for the two-loop self-energy correction for the 
ground state of hydrogen-like ions with with the nuclear 
charge number Z = 10 — 30. The higher-order (in Za) 
remainder function is inferred from the numerical results 
and extrapolated towards Z = and 1. The extrapolated 
value is in marginal agreement with the analytical result 
obtained within the perturbative approach. 

The work presented in this paper was supported by 
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Appendix A: Dirac Coulomb Green function in the coordinate-momentum representation 



The Dirac Coulomb Green function is commonly written in coordinate space as an expansion in the relativistic 
angular momentum parameter k [11, 12, 29], 

r(F -r -r \ = V ( G « 1 ( jB ' a;i ' a; 2)XK ([1 (*i)x^(*2) Gl 2 (E, Xl ,x 2 ) {-i)x K »(xi)X- K »(x2) \ fAn 



where XKfi (x) are the Dirac spin-angular spinors [23] and 
x = x/ 1 x I . The 2x2 matrix of the radial components G l J 
is referred to as the radial Green function and denoted 
as G K . The radial Green function can be expressed in 
terms of the two-component solutions of the radial Dirac 
equation regular at the origin ((/>°) and the infinity (<^£°) 
as follows 

G K {E,x u x 2 ) = - ( j)™{E,x 1 )4>f{E,x 2 )6{x 1 -x 2 ) 



(E,x 2 )9(x 2 -xi). 

(A2) 



The upper and the lower components of the functions 



(jP K and </>£° will be denoted by subscripts "+" and 
respectively. They are given by 



x Z/2 

aZ 



(A-i/)M„ i x (2cx) 



— ) M v+hx (2cx) , (A3) 



c ± (m=a-v 2 ^M 

±W u+lx (2cx) 



(A4) 
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where A K = 4c 2 T(l + 2X)/T(X-v), c = y/l - E 2 defined 
so that Sft(c) > 0, A = ' n 2 - (Za) 2 , v = ZaE/c, and 
M ai p and W a ^ are the Whittaker functions of the first 
and the second kind, respectively. 

The Dirac Coulomb Green function in the coordinate- 
momentum representation is obtained from the above for- 
mulas by the Fourier transform over one of the radial ar- 



J 



guments. Let us consider the transform over, e.g., the 
second radial argument, 



G(E,X! 



■*■>-/ 



dx 2 e ip2 - X2 G(E, Xl ,x 2 ). 



Its partial-wave expansion takes the form 



(A5) 



(A6) 



C'(E x v ) = V?* [ G « 1 (- B . a; i.P2)x«M( i i)xl eM (P2) G 12 (E,x 1 ,p 2 ) X n,,(x 1 )x- KI ,(P2) \ 
' 1 ir \Gl 1 (E,x 1 ,p 2 )ix- K ^i)xUp2) Gf{E 1 x llP2 ) lx ^{x 1 ) x l K ,(P2) ) ' 

with the radial part given by the following matrix 

r .„ , . [°°, 2 ( 3i(P2X 2 )G?(E,x 1 ,x 2 ) -^j T (p2X 2 )Gi 2 (E,x u x 2 )\ 
G "i E >*>P>)= 4 *J o dx 2 x 2 y jl{p2X2)Gf{E , XuX2) - jj(p 2 x 2 ) Gl?(E, xi,x 2 ) ) ' 

where I = \n + 1/2| — 1/2 and 1 = \k — 1/2| — 1/2. Using Eq. (A2), we obtain the following representation for the 
radial Green function in the mixed coordinate-momentum r epresentation, 



(A7) 



G K (E,xi,p 2 ) 



-4%{E,xi)1>f{E,p2;xx) 
$(£,zi)^~ T (E,P2;a;i), (A8) 



where 



ip°(E,p; xi) =4tt / dx 2 x\ 



and 



^~(S,p;a;i)=47r 



3i{px2)4>° K „+{E,x 2 ) 
-]%[jj(px2)<t>° K ,-(E,x 2 ) 



/ C?X2 ^2 

I Xi 

^-^.n( P X2)^-(E,x 2 ) 



(A9) 



. (A10) 



The integration over x 2 in the functions ip® & n d 
has to be performed numerically. The problems here are 
that (i) the integration interval depends on x\ and (ii) the 
integrand contains the spherical Bessel function which 
oscillates rapidly in the high-momenta region. Clearly, a 
straightforward use of Eqs. (A9) and (A10) in our cal- 
culations would lead to a re-evaluation of the integral 
for each new value of X\, making the computation pro- 
hibitively expensive. One can observe, however, that if 
the function ip K (E, p; x) is known for a particular set of E, 
p, and x, then the evaluation of ip K (E,p;x') can be done 
by computing the Bessel transform integral over the in- 
terval (x, x') only. So, introducing an ordered radial grid 
{x^, one can store the whole set of values {ip K (E,p; Xi)} 
by performing just one Bessel transform over the interval 



(0,oo). This shows that for a fixed values of E and p, 
the integrations of the type L dx f(x)ip K (E,p;x) can 
be performed without a recalculation of the Bessel trans- 
form integral. 



Appendix B: Free Dirac Green function in the 
coordinate-momentum representation 

The free Dirac Green function G^ is a much simpler 
object than the Dirac Coulomb Green function G and 
is known in the closed analytical form as well as in the 
partial-wave expansion form (see, e.g., Ref. [12]). For 
the purposes of the present investigation, we employ the 
coordinate-momentum representation and put G^ into 
the form analogous to Eqs. (A6) and (A7). The simplest 
way to achieve this is to start with the momentum repre- 
sentation of G(°) , which has a particularly simple form, 



G^(E, Pl ,p 2 ) =(2tt) 



3 5 i (p 1 -p 2 ) 
j°E — 7 • p 2 — m 



E + a ■ p 2 + m7° , 
( 27r ) ^ :l „2 * V " Pa) > 



E 2 - p\ - m 2 



(Bl) 



where a = 7%. Using the completeness of the angular- 
momentum spinors Xk^i, 



2X"m(Pi)x^(P2) = IS(p!-p 2 ), 



(B2) 



and the identity (er • p)x^{p) = -X-k»{p) , we cast 
Eq. (Bl) into the partial- wave expansion form similar to 
that for the Dirac Coulomb Green function, 
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E 2 -p 2 2 -m 2 ^y -P2X-^{Pi)xl„(P2) {E-m) X -^{pi)X-^{P2) 



(B3) 



where Pi = \pi\. The coordinate- momentum representation of is obtained by the Fourier transform of the above 
expression over the first radial argument, 

G^(E,x u p 2 ) = J -^Le^G (0 HE, Pl ,p 2 ). (B4) 

After performing the integration over pi, the free Dirac Green function is written in the form of Eq. (A6), with the 
radial part given by 

r(o) ( j? ^ „\ 4?r ( (E + m)ji(p2Xi) ~P2ji(P2Xi)\ . . 

G « {E ^ P2) - E^-pl-m 2 { fttoifax!) -^(E-m)jj(p 2Xl ) ) ■ 



Appendix C: Angular factors 



and r k 



In this section we address the factors t Kri . K , . , . 

which are defined by Eqs. (35)-(37). Inserting the explicit 
definitions of the angular-momentum spinors in these for- 
mulas, averaging over the momentum projections of the 
reference state, and calculating the sums of the Clebsch- 
Gordan coefficients, we arrive at the following results, 



^Kn ,K a (<-0 



(_l)J.+l/2 n . n 



47T 



IE 



ja jn J 



'Z(-l) M Yj M (q)Y£ l - M (p 1 ,p 2 ), (CI) 



M 



(_l)jq + l/2 U Jn 



M V .JM 



J 

1/2 

■J -Mi 



Y^(q,P2)YS nl -™(p u p 2 ), 



M 



(C4) 



/ ft IT. „ I Ja J" J 

Cko (JL) = (-1)^JA^£ 1/2 1/2 1 

V ^ U i« L (la In L 

x ^(-l) M Y LM (q)Y l L nl - a M ( Pl ,p 2 ), 

M 



(^) = 



(-i) J ° +1/2 n jn f Jo i„ j 



(C2) 



V3 % \ i„ 1/2 

(C3) 



A/ 



where 

^ J , M (pi,p 2 ) ' 



are the 



v /(2.n + l)(2j 2 + l)... and 
bipolar spherical harmonics 



iil 2 
[25]. 

With help of formulas from the book [25], it is possible 
to obtain explicit results for the angular factors £ K „ jKa 
and s£ Ka , which are functions of pi = [pi \ , p 2 = \p 2 | , 
and q = \pi — p 2 \ only. However, the resulting formu- 
las turn out to be rather lengthy and not very conve- 
nient for numerical evaluation as they become numeri- 
cally unstable for q — > 0. Because of this, we prefer to 
evaluate Eqs. (C1)-(C4) numerically, after some simpli- 
fications that exploit the fact that the result does not 
depend on any angles except for pi ■ p 2 . Namely, we set 
the azimuthal spherical coordinate of pi and p 2 to zero 
(0i = 4> 2 = 0) and direct pi along the z axis (6\ = 0). 



